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Abstract 

This paper addresses the robust counterparts of optimization problems containing sums of 
maxima of linear functions. These problems include many practical problems, e.g. problems with 
sums of absolute values, and arise when taking the robust counterpart of a linear inequality that 
is affine in the decision variables, affine in a parameter with box uncertainty, and affine in a 
parameter with general uncertainty. 

In the literature, often the reformulation is used that is exact when there is no uncertainty. 
However, in robust optimization this reformulation gives an inferior solution and provides a pes¬ 
simistic view. We observe that in many papers this conservatism is not mentioned. Some papers 
have recognized this problem, but existing solutions are either conservative or their performance 
for different uncertainty regions is not known, a comparison between them is not available, and 
they are restricted to specific problems. We describe techniques for general problems and com¬ 
pare them with numerical examples in inventory management, regression and brachytherapy. 
Based on these examples, we give recommendations for reducing the conservatism. 


1 Introduction 

Robust Optimization (RO) first appeared in Soyster (1973), and after receiving very little attention 
in the subsequent decades it has been an active research area since Ben-Tal and Nemirovski (1999) 
and El Ghaoui and Lebret (1997) started publishing new results in the late nineties. In its basic 
form, it requires a solution of an optimization problem to be feasible for any realization of the 
uncertain parameters in a given uncertainty set. For several choices of the uncertainty region this 
leads to tractable problems; for instance the robust counterpart (RC) of an LP with polyhedral 
uncertainty can be reformulated as an LP and the RC of an LP with ellipsoidal uncertainty is a 
conic quadratic program (CQP) (Ben-Tal et ah, 2009a, p. 21). 

In RO two distinct formulations that are equivalent in the nonrobust case may have different 
RCs. This is the case for optimization problems containing the sum of maxima of linear functions, 
of which the sum of absolute values is a special case (|x| = max{x, —x}). These problems arise in 
inventory management, supply chain management, regression models, tumor treatment, and many 
other practical situations. 
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Because RO is applied constraint-wise w.l.o.g. and the objective function can always be formu¬ 
lated as a constraint (Ben-Tal et ah, 2009a, p. 11), we focus on the following robust constraint: 

max{£ij(C, x)} < d \/C e Z, (1) 

i^i 

where £ and lij are biaffine functions in the uncertain parameter G and the decision variable 
X G M”, d G M is the right hand side, and Z C is a user-specified uncertainty region, e.g. a box 
or an ellipsoid. We can assume this uncertainty region to be closed and convex w.l.o.g., as the left 
hand side is convex in ^ and consequently the worst case is always located at an extreme point of 
the uncertainty region. In the remainder of this paper our objective is always to minimize d, but 
all methods can still be applied when a different objective is used. 

In the literature, often the following RC is used, obtained from a reformulation with analysis 
variables y* that is exact when there is no uncertainty: 


(RC-R) t{Cx)+Y,yi<d 

VC G ^ 

(2) 

iei 



Vi > ^ij iC,x) 

VC G ^ Vi G I Vj G J. 

(3) 


In this formulation, to which we refer as the RC-R (RC of the reformulation), y* G M is a fixed 
variable, taking the worst case value of the term of the sum. In many cases the terms of the 
sum do not all reach their worst case value in the same realization of the uncertain parameter (), 
and therefore the RC-R is a conservative reformulation of constraint (1), i.e. a solution (x,d) that 
is feasible for RC-R is also feasible for constraint (1), but not necessarily vice versa. However, this 
reformulation is frequently used without mentioning its conservatism. It has the advantage that 
the constraints are linear, so that tractable reformulations exist for many uncertainty regions. 

Bertsimas and Thiele (2006) use the RC-R for a robust inventory problem which includes the 
sum of holding and backlogging costs: max{cftXj, —CbXj}, where Xi is the inventory level at 

time period i and Ch and are the holding and backlogging costs per time period respectively. 
Their uncertainty region is the intersection of {(" : ||Ci:i|li < Bj} Vi and {() : HClIoo — where 
Ci:i is the vector consisting of the first i elements of C, and Tj is not necessarily integer (budgeted 
uncertainty). Their formulation allows that the values of different analysis variables take values 
based on different realizations of the uncertain parameter. The same uncertainty region is also 
used by Bertsimas and Thiele (2004) for a supply chain model, by Thiele (2004) for supply chains 
and revenue management, by Alem and Morabito (2012) for a production planning problem under 
uncertain demand, and by Wei et al. (2009) for a slightly more complicated inventory model in 
which items can be returned and remanufactured or disposed. With the exception of (Thiele, 2004, 
p.37), none of these papers mentions the conservatism of the formulation. 

Ng et al. (2010) treat a lot allocation problem where each order is assigned to one or more pro¬ 
duction locations before the production capacity of each location is fully known. Every production 
location is assigned to at most one order. Their uncertainty region is an ellipsoid over all production 
locations, while an analysis variable is used for every order, so their formulation is conservative. 

Kropat and Weber (2008) consider a robust linear cluster regression model using the sum of 
absolute values with polyhedral and ellipsoidal uncertainty regions. They introduce an analysis 
variable for every absolute value, which is conservative because every uncertain parameter affects 
two absolute values. 

Ben-Tal et al. (2005) solve a multi-period inventory problem and use affine decision rules for the 
actual decisions (AARC). For every time period there are analysis variables indicating the costs in 
that time period. These costs are given by the maximum of holding costs and shortage costs, which 
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both are linear functions of past demand. They note that an analysis variable should therefore be 
replaced with a function being the maximum of the two linear functions. Because that would lead 
to a very complicated robust counterpart, they replace the analysis variables with linear decision 
rules instead, which is conservative. Replacing analysis variables with linear decision rules was also 
done by Ben-Tal et al. (2009b) and Ben-Tal et al. (2011). 

Bienstock and Ozbay (2008) were the first to identify and eliminate the conservatism of the 
RC-R. The idea behind their solution is that it suffices to make constraint (1) hold for just the 
vertices of the uncertainty region. The constraint then also holds for all other elements in the 
(convex) uncertainty region, because the constraint is convex in C. They generalize it to a cutting 
plane method that only adds a subset of the vertices, which they successfully tested for computing 
basestock levels under budgeted uncertainty. It is not yet known how well this cutting plane method 
performs on other problems or different uncertainty regions. 

In this paper, we make the following contributions. First, we identify that a conservative 
reformulation is often used in the literature without mentioning its conservatism. Second, we make 
a classification of solution approaches. Third, we present a new cutting plane method (Algorithm 
2). Fourth, we show that a linear constraint with biaffine uncertainty is a special case. This 
provides a method to solve or approximate a conic quadratic constraint with box uncertainty, for 
which currently no efficient methods are available. Fifth, we demonstrate all approaches on several 
numerical examples. Such a comparison is currently not available. Sixth, based on our numerical 
examples, we give recommendations for reformulating inequalities containing sums of maxima of 
linear functions. 

The structure of this paper is as follows. Since it is not possible to compare the objective value 
of different formulations in RO, we introduce a new performance number which is independent of 
the reformulation in Section 2. In the examples we consider, this performance number is the slack in 
constraint (1). We provide an overview of exact (non-conservative) formulations, approximations, 
and cutting plane methods in Section 3. In that section we show similarities between different 
methods, and also show how several approaches can be combined. The application scope of this 
paper is extended in Section 4, where we show that our methods can be applied to an uncertain 
conic quadratic constraint with box uncertainty. We evaluate the methods in Section 5 on some 
small toy problems and three larger problems. We give conclusions in Section 6, and show how the 
robust counterparts in the aforementioned papers could be improved. 


2 The true robust value 


This section explains how two solutions from different reformulations can be compared. For a 
minimization problem containing sums of maxima of linear functions, we define the true robust 
value of a solution to be the maximum over (( in ^ of the unreformulated problem with all decision 
variables fixed. We have assumed that in this paper our objective is always to minimize d, an 
analysis variable at the right hand side of (1), so in our case the true robust value of a solution x 
is: 


Vtrue{x) = max < £(C,x' 


-h ^max{£y(C,x)} 
is/ 


(4) 


Determining this value is a difficult problem, because it requires the maximization of a convex 
function over a convex set. The global maximum over ^ G ^ of the sum of maxima of linear 
functions does not necessarily coincide with a maximum over G ^ of one of those linear functions. 
One way to obtain the exact value is by solving the following optimization problem with integer 
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variables for fixed x: 


max 

s.t. 


+ (5) 

i£l 

Vi < ^ijic, x) + M(1 - Zij) yiG I yj e J 

Zij = 1 yi G I 

j&J 

CGZ,yG e 


where M is a sufficiently large number. This problem is an MILP for polyhedral uncertainty, and 
an MIQCP for ellipsoidal uncertainty. 

Another way of obtaining the exact value is by considering linear optimization problems, 
e.g. for I = {1,2} and J = {1,2} consider the following four optimization problems: 

max£(C, x) + 4,i(C, x) + £2,i(C, x), (6) 

max£(C, x) + x) + £2,2(C, x), 

inax£((,x) +£i, 2 ((,x) + £ 2 ,i(C,x), and 

max£(C, x) + 4,2(C, x) + £2,2(C, x). 

Each of these problems is easy, because the maximum of an affine function over a box or an ellipsoid 
can be computed in a few operations. The true robust value is the largest value of the computed 
maxima. 

If determining vtrue(x) is intractable, bounds can still be obtained. Filling in any ^ from the set 
Z at the right hand side of equation (4), for instance the nominal value, gives a lower bound. If the 
dimension of ^ is small, the bound can be improved with global optimization techniques. Upper 
bounds can be obtained by fixing x in any conservative reformulation (such as those mentioned in 
Section 3.2) of constraint (1) and optimizing over the other variables. 


3 Solution approaches 

This section lists exact solution approaches and approximation methods, most of which can be 
applied to general RO problems containing the sum of maxima of linear functions. Many of these 
methods have been used before, but this full classification is new. This allows us to show the 
similary between some methods, and to show how approximations can be combined with exact 
methods. 

3.1 Exact reformulations 

3.1.1 Vertex enumeration 

Vertex enumeration is an exact solution method first used for a RO problem containing the sum 
of maxima by Bienstock and Ozbay (2008) that is powerful especially when the uncertainty region 
has a small number of vertices. Let V denote the finite set of vertices, and consider the following 
reformulation of constraint (1): 


(Vertex enumeration) £((, x) + yj<d 

VC G u 

( 7 ) 

i£l 



yi > £ij{C^x) 

yiGi yj gj VC g V. 

(8) 
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Constraint (8) is no longer a semi-infinite constraint, because C is a finite set. This reformulation 
is exact because the left hand side of constraint (1) is convex in and a convex function takes its 
maximum at an extreme point of its domain. 

3.1.2 Enumeration of robust linear constraints 

The RC-R is inexact because it has analysis variables that may take values corresponding to different 
worst case scenarios. A constraint with a single max function does not suffer from this problem, 
because an equivalent set of linear constraints can be formulated without analysis variables. An 
exact reformulation of RC (1) can be obtained by first rewriting it as a constraint with a single 
max{-} function by enumerating all combinations, and then applying RO to the reformulation: 

(EORLC) f(C,x) + ^%q(C,x) <d Vj(i) G J VC G 

i&l 

We call this the enumeration of robust linear constraints (EORLC) formulation. It has the ad¬ 
vantage that the constraints are linear, so that tractable reformulations exist for many uncertainty 
regions. For example, with I = {1,2} and J = {1,2}, the EORLC formulation has the following 
constraints: 


f(C, x) -h 4,i(C, x) + £2 ,i(C,x) < d VC e z 

£(C,x) -h 4,i(C, x) + £2,2(C,x) < d VC e z 

£(C, x) + £i,2(C, x) + ^2,i(C, x) < d VC e Z 

£(C, x) + £i, 2 (C, x) + £ 2 , 2 (C, x) < d VC € z. 


A similar formulation is also given by Bienstock and Ozbay (2008), where it was neglected for 
its exponential size. While the number of constraints | J|VI indeed grows exponentially with the 
number of terms in the summation, it is effective for small |/|. There are situations in which I is 
small indeed, for instance with a planning horizon of up to ten periods. 

EORLC has a strong relation to vertex enumeration that was not observed before to the best 
of our knowledge. In fact, EORLC is vertex enumeration on a different set. Constraint (1) can be 
formulated as: 


VAi € AVI-1 VC G 

ieijeJ 

where = {Aj G M^l ; J2j^j ^ij = standard simplex in The vertices of 

the simplex are given by unit vectors. If A* is a unit vector, JCjej ^ij£ij{Cjx) simplifies to a single 
£ij{C,x). It follows that vertex enumeration on the |/| simplices gives the EORLC reformulation. 

EORLC can benefit from two preprocessing steps in order to reduce the final number of con¬ 
straints. First, every max{-} term should contain at most one function that does not depend on C- 
If it has more than one, those functions can all be replaced with a single analysis variable. Second, 
inequalities of the form amax{0,£{C,x)} + 6max{0, —£{C-,x)} < d with a, 6 > 0, which are often 
used for holding and backlogging costs, should be reformulated as max{a£(C,x), —b£{C-,x)} < d. 

3.1.3 Cases with special structure 

There are several special cases of (1) that allow an exact reformulation. This section lists a few 
general cases. More specific cases can be found in (Ben-Tal et ah, 2009a, Ch. 12.2) and in Xu et al. 
(2009). 
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The first case is when the uncertainty region Z is the direct product of sets Zi^ where term i in 
the left hand side of constraint (1) is only affected by Zi. The RC-R is exact because all analysis 
variables can take their worst case values simultaneously. 

The second case is when the inequality contains the sum of absolute values of linear functions 
of Ci and x\ 


^{C,x) + '^\ai{x) +Pi{x)\\<d yC^Z, (9) 

i£l 

where a* : MZ —)■ M and /3i : —)■ M are linear functions, the components of ,0* that may be nonzero 

for one i are zero for all other i in I, and the uncertainty region Z is centrosymmetric around C = 0, 
i.e. Z is closed under changing the sign of one or more vector elements (a box and an ellipsoid are 
examples of such sets). Note that the assumption that the symmetry is around 0 is made w.l.o.g. 
The following constraint is equivalent to (9): 

^{C,x)+^{\ai{x)\ +Pi{x)Ci} < d yC^Z. (10) 

i£l 

Equivalence is readily checked by conditioning on the sign of ai{x). The formulation for (10) where 
absolute values are replaced with analysis variables is equivalent to (10), since the analysis variables 
do not depend on (. 

The third case is when, for a fixed i £ I, each of linear functions under the maximum is the 
sum of a common nonnegative linear function of C and a linear function of x: 

^(C, 2 :) + 51 niax{ai(C) + /3i{()iij{x)} <d VC G 2:, 
ie/ 

where Oi : —)■ M and j3i : —)■ M+ are linear functions. The common functions of C can be 

placed outside the max{-} expression: 


£(C,x) + 


iei 


oiiiO 


A(C) max{£ij 

j&J 


(^)} 


< d 


yCGZ, 


and the RC-R of this constraint is exact. If the range of jdi is M instead of M+, then maxjgj{£jj(x)} 
should be m.mj^j{iij{x)} when /0i(C) < 0. This can be modeled as follows. Given a subset /+ C /, 
we define a set which consists of those C for which ,0i(C) > 0 for i £ and /3i(C) < 0 for i £ /\/+: 


Z(/+) = zn{C:/3i(C) >0 VfG/+, /3*(C)<0 yi£l\I+}. 

Note that 

Z= U Z(I+). 

Rc/ 

The constraint can now be written as: 

/^i(C)^ax{£ij(x)}+ /^i(C)inm{£ij(x)} < d yC £ Z(I+) V/+ C /. 

i£l ie/+ 

The number of constraints is 2^1 (one for each /+ C /), which is less than the constraints 
obtained with EORLC. The max and min expressions do not depend on C, so the reformulation 
with analysis variables is exact. Each constraint is still convex despite the min expressions, because 
their coefficients /di{() are negative. 
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3.2 Conservative approximations 

The RC-R (2)-(3) is a conservative approximation to (1). We discuss how the conservatism can be 
decreased, and also present a new method. 

The variables y* in the RC-R have been introduced for modeling the max{-} expression. We do 
not need to know their values because they do not correspond to a “here and now” decision. Only 
X has to be known for implementing a solution. The values of yi may be adjusted according to the 
realization of the uncertain parameter ^ as long as the constraints hold for every realization of C 
in the perturbation set (Adjustable RC). As first applied by Ben-Tal et al. (2005) for a multi-stage 
problem and later also done by Ben-Tal et al. (2009b) and Ben-Tal et al. (2011), we can make yi 
an affine function of C, which leads to the Affinely Adjustable RC of the reformulation (AARC-R). 
After substituting yi = Vi + (with decision variables Vi € R and Wi G M^), constraints (2)-(3) 
become: 


(AARC-R) £(C, a:) + ^ (ui + rcjc) < d VC G ^ 

i&I 

Vi + wlc > x) yc e z \/i e I \/j g J. 

This substitution gives a less conservative reformulation, while the robust counterpart is often still 
tractable, because robust linear constraints are tractable for a wide class of uncertainty regions. 
The power of the AARC-R can be increased by lifting the uncertainty region to a higher dimension 
(Chen and Zhang, 2009). If ^ij does not depend on one or more components of C for all j for some 
fixed i, the computational complexity can seemingly be reduced by making yi a function of only 
those components of C that appear in term i for one or more j. However, it is easy to construct an 
example where this reduction introduces more conservatism. 

There are two different approaches that also lead to the formulation AARC-R. The first approach 
is considering the Fenchel dual problem of maximizing the left hand side of constraint (1) over ^ 
in Z. We give the full derivation and a proof of equivalence to the AARC-R in Appendix A. The 
other approach is derived in Appendix B. 

An approach that is less conservative than an affine decision rule, is a quadratic decision rule: 

yi = Vi + w]C + CWiC, 

where Vi G M, rcj G and Wi G R^^^ are new analysis variables. This is called a Quadratically 
Adjustable RC of the reformulation (QARC-R) which is known to be tractable for ellipsoidal 
uncertainty and (under some restrictions) for box uncertainty. A deterministic reformulation of the 
QARC-R with ellipsoidal uncertainty is given in Appendix C, resulting in |/||J| -|- 1 LMIs of size 
IL -|- 11. For box uncertainty, when the quadratic terms are restricted such that each element of C is 
multiplied with itself and at most one other element of Q we can use the result by Yamkoglu et al. 
(2012) to write the RC as an SDP with (|/||J| -|- l)|'T/2] variable matrices of size three. 

A new way to reduce the conservatism of the RC-R is by first combining several max expressions 
before reformulating. In order to do this, partition I into |G| groups: I = U^eG where the groups 
are mutually disjoint: Ig^ n 1^2 = 0 51 > 1/2 £ G,gi / 52 > and partition in such a way that \Ig\ is 

small for all g. Now introduce analysis variables for all g: 

^{C,,x) + ^yg<d yC,eZ 

g&G 

yg>Yl niax{Qj(C,x)} e Z yg eG. 

i&Ig 
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Because the cardinality of Ig is small, the sum of max expression within each constraint can be 
transformed into a single max (EORLC). Each constraint in this reformulation is therefore tractable 
and can be solved exactly, but conservatism still comes from the analysis variables yg. These analysis 
variables can also be written as a linear, quadratic or more general function of C- 

3.3 Cutting plane methods 

In this section we describe two cutting plane methods. The first one is based on vertex enumeration 
and was used in Bienstock and Ozbay (2008), the other one is new and based on EORLC. 

Vertex enumeration results in very large problems if there are many vertices. Eor box uncer¬ 
tainty, the number of vertices grows exponentially in the time horizon. Also for budgeted uncer¬ 
tainty, which is described in Section 1, the number of vertices quickly becomes very large. The 
cutting plane method outlined in Algorithm 1 adds only a subset of the vertices. 


Algorithm 1 Cutting plane method based on vertex enumeration 
Require: A linear program LP with constraints (7) 

1: V := (the nominal value) 

2: k :=0 
3: repeat 
4: k := k + 1 

5: Solve LP 

6: Let X* be the minimizer of LP, and LB be the value of d in LP 

7: Let UB := vtrue{x*), and be its maximizer 

8: V:=VU{C^} 

9: until UB — LB < e 


While the algorithm is running, LB is a lower bound on the optimal value of d in (7) because 
it is the value of a relaxation of the original problem, and UB is an upper bound on the optimal d, 
because it is the maximum for some x and not for the optimal x. The difference UB — LB indicates 
the current violation of the constraint. If e is set to a larger value, the algorithm terminates 
quicker but does not give the optimal solution. If e = 0 and the algorithm terminates, the final 
solution is robust feasible and robust optimal. It is possible that this algorithm still enumerates all 
constraints, but we have not encountered problems for which this is the case. Bienstock and Ozbay 
(2008) find that for their problem the number of iterations does not increase when the time horizon 
is increased, and is around four on average even for 150 time periods. The number of vertices of 
their uncertainty region cannot be deducted from their report. In our numerical examples we often 
find a larger number of iterations. The same cutting plane method was used by Bohle et al. (2010), 
again for budgeted uncertainty, but solution times and the number of generated constraints are 
not reported. It is still an open question how well this method works on other problems, and for 
non-polyhedral uncertainty regions. 

The most time consuming step is determining the true value, since this involves the maximiza¬ 
tion of a convex function. As pointed out by Bienstock and Ozbay (2008), it is not necessary to 
find the optimal solution. It suffices to find any C, for which the maximization problem has a larger 
objective value than LB, because it corresponds to a violated constraint in the relaxation, and 
adding that constraint strengthens the relaxation. Bienstock et al. report that also the problem 
LP does not need to be solved to optimality. Optimization can stop as soon as the objective value 
is less than UB, because at that point the current solution is already better than the solution in 
the previous run. 





The second cutting plane method is new and based on EORLC. If |/| becomes too large for 
EORLC, this cutting plane method adds only a subset of the robust linear constraints. It starts 
with the nominal problem and adds new robust constraints until the solution is robust feasible 
(Algorithm 2). 


Algorithm 2 Cutting plane method based on EORLC 
Require: A linear program LP with constraints (7) 

1: V := (the nominal value) 

2: k :=0 
3: repeat 
4: k := k + 1 

5: Solve LP 

6: Let X* be the minimizer of LP, and LB be the value of d in LP 

7: Let UB := vtrue{x*), and be its maximizer 

8: Let jf be the maximizer of max^gj x) 

9: Add the following robust constraint to LP: 

^(C, 4^) + 51 (‘5^’ x)<d z 

iei 

10: until LfB — LB < e 


Similar to Algorithm 1, a lower and upper bound are reported and the stopping criterion can 
be adjusted. Also for this algorithm, we have not encountered numerical examples in which all 
constraints are enumerated. The main difference between this algorithm and Algorithm 1 is the 
constraint that is added in every iteration. 

Just as with the cutting plane method based on vertex enumeration, determining the true value 
is the most time consuming step, and also here the optimization problem for determining the true 
value does not need to be solved to optimality as long as the value is larger than LB. Also LP 
does not need to be solved to optimality. 

Both cutting plane methods can easily be combined so that two sets of constraints are added 
per iteration. This is done by inserting line 8 of Algorithm 1 in between lines 9 and 10 of Algorithm 

2 . 

The stopping criterion in both algorithms is UB — LB < e. It is also possible to use a dimen¬ 
sionless stopping criterion such as 2{UB — LB)/{I + \UB + LB\) < e. This stopping criterion is 
based on the relative gap. 


4 RC of a linear constraint with biaffine uncertainty 


Constraint (1) may appear in RO itself when a constraint has biaffine uncertainty, and the uncer¬ 
tainty region of one parameter is a box. In general, such a constraint can be written as: 










( 11 ) 


where I: xM^ —)■ M is a triaffine function, p is the radius of the box, and Z is the uncertainty 

region of Constraints with biaffine uncertainty have never been investigated before to the best 
of our knowledge. To show the equivalence to constraint (1), choose i, ii and I in such a way that 


9 





constraint (11) is equivalent to: 

<(C®,:=) + EC“’<.(C®.^)<<i vc<»: C« <p vc^ez, 

i^I 

and maximize the left hand side with respect to 

l{C,^^\x)+Y^p\l,{C^'^\x)\<d VC(2) GZ, (12) 

iei 

which is indeed in the form of constraint (1). This RC includes many practical problems, of which 
we give four examples. 

The first example is a constraint a x < d where each element of a has been estimated using 
a regression model with the same regressors Q for every element: Oj = ,0? + (/3^*^)^C + {Pf ^ 
G R^,ei G M), in that case the constraint is: 

+ BC + syx < d. 

We assume that the coefficients G and B G the matrix having as rows, have been 

estimated with some error, and also that the regressors C are not fully known. The constraint can 
be written as constraint (11) if /3®, B, (( and e lie in some specified uncertainty region, and the 
uncertainty of either R or is a box. 

The second example appears in Adjustable RCs. Consider a robust constraint: 

^{y^\y‘^\x) + b’y <d <p e Z, 

OO 

where y G is an adjustable variable that represents a “wait and see” decision that can be made 
after and are (partially) observed, and 6 is a vector of known coefficients (fixed recourse). 
Determining the true optimal policy for y as a function of and is often intractable, which 
is why a suboptimal y is often determined by a parameterization, i.e. y is written as a function 
of which the coefficients are decision variables. The first parameterization in RO was proposed by 
Ben-Tal et al. (2004) who proposed an affine decision rule, which results in a problem which is in 
the same class (LP, CQP or SDP, depending on the uncertainty region) as the problem with y as 
a “here and now” decision variable. This was later extended to a quadratic decision rule with an 
ellipsoidal uncertainty region, resulting in an SDP (Ben-Tal et al. (2009a)), and to a polynomial 
decision rule of arbitrarily large degree restricted to uncertainty regions described by polynomial 
inequalities, resulting in a conservative SDP (Bertsimas et al. (2012)). The latter includes the 
biaffine decision rule y = i where £ is a function x x _s, and v is 

a vector of coefficients to be determined by the model. This decision rule could be very useful if 
a problem is affected by several sources of uncertainty. The advantage over affine decision rules is 
that the former includes cross terms of and Applying the results of Bertsimas et al. (2012) 
gives an SDP which is not only conservative, but also has a potentially large instance size. E.g. if 
|I| (the dimension of C^^^) is very small but is in a box of (large) dimension L, the conservative 
SDP has at least one variable matrix of size (|/| -|- 1)T-|- |/| -|-1 and over 2|L| smaller matrices, while 
our result is exact and gives a practically solvable LP. If instead is in an ellipsoid of (large) 
dimension L, the SDP still has the large variable matrix of size (|/| -|- 1)L -|- |/| -|- 1 while our exact 
result gives a CQP for which efficient solvers are available. 

The third example is a constraint with unknown coefficients and implementation error. Consider 
the following constraint: 

< d g ^i, (13) 
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where £ is a vector of n functions that are linear in the uncertain parameter Now suppose that 
there is implementation error in x, i.e. instead of x we implement a vector of which component i 
is given by Xi + Q ' (additive implementation error) or by Q Xi (multiplicative implementation 
error). After substituting x, constraint (13) becomes: 

£(C^^^)"(x + C^^^) < d G G ^2, 

in case of additive implementation error, and: 

n 

< d G G ^2, 

i=l 

in case of multiplicative implementation error. Both constraints are special cases of constraint (11) 
if either Zi or Z 2 is a box. 

The fourth example is the following robust constraint with box uncertainty: 

^(C^^\a:) + XI < d < P, (14) 

where is a vector of L linear functions, and I equals 1, 2 or 00 . In order to see that this constraint is 
equivalent to constraint (11), note that constraint (14) is a reformulation of the following constraint: 

£{c^^\x) + Y.{ci^^yh{&\x)<d <p cf *<i ykeK, 

k^K 

where H-H^ is the dual norm, and is a vector in for all k in K. We will show how the results 
in this paper can be used for solving the robust constraint (14) for different choices of 1. 

For I = 2 and |Ar| = 1, constraint (14) is a robust conic quadratic constraint, for which the 
RC is known only in special cases, one of which is when the vertices of the uncertainty region can 
be enumerated (Ben-Tal et ah, 2009a, p. 159). The case |ih > 1 has not been covered yet. A 
(conservative) reformulation with analysis variables for every ik{y^\x) reduces the problem to 
a problem with one linear constraint and \K\ robust conic quadratic constraints, each of which is 
of the form (14) with \K\ = 1, so it can be reformulated using vertex enumeration. A different 
approach, that is not only exact but also results in a smaller problem than obtained by using analysis 
variables, is to apply vertex enumeration to constraint (14) directly. Vertex enumeration is exact 
because the constraint is convex in If the dimension of the box is very large, (iterative) vertex 
enumeration is no longer tractable, so even the case |Ar| = 1 becomes unsolvable, and tractable 
conservative reformulations are not known. Our reformulation (12) allows the use of all approaches 
from Section 3. 

For I = 1 constraint (14) is a special case of constraint (1), and we know from Section 3 how to 
solve it or how to find a conservative reformulation. The reformulation may be useful if it is solved 
with the RC-R, the AARC-R or the QARC-R, because the resulting formulation is very different. 
Vertex enumeration and EORLC are not useful on the reformulation (12), because they correspond 
with EORLC and vertex enumeration on the original constraint, respectively. 

Eor I = 00 constraint (14) is a special case of constraint (1) with I = K and J = {1,2, ..,T}, 
and we know from Section 3 how to solve it or how to find a conservative reformulation. The 
reformulation (12) may be useful if it solved with the RC-R, the AARC-R or the QARC-R, because 
the resulting formulation is very different. The reformulation allows to do vertex enumeration on 
which may be faster than vertex enumeration in the original constraint (which is done on the 
vertices of C^^^) if L and |Ar| are small relative to |/|. 
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For all three cases of /, it holds that when EORLC is used on the reformulation (12), the 
resulting constraints are the same ones resulting from vertex enumeration on constraint (14). This 
implies that reformulating (12) is a redundant step if it is followed by EORLC. 

5 Numerical examples 

The LP, MILP, CQP and MICQP problems in this paper have been solved with AIMMS 3.11 PR3 
x32 with ILOG CPLEX 12.1 unless stated otherwise. SDP problems have been modeled with CVX 
(Grant and Boyd (2010)) and solved with SDPT3 (YALMIP (Lofberg (2012)) would nowadays be 
a better choice as it allows the user to specify a linear constraint with quadratic uncertainty with 
an ellipsoidal uncertainty region directly). Reported computing times have been obtained under 
Windows XP SP3 on an Intel Core2 Duo E6400 (2.13 GHz) processor and 2 GB of RAM. 

5.1 Computing the true robust value 

We have run some experiments to determine how quickly the true robust value can be computed 
using the two different methods from Section 2. Using the problem with integer variables (5), the 
worst case C can be determined in less than a minute for |/| =50 and | J| = 3 for box uncertainty 
of dimension 50, and in 7 seconds for |/| = 20 and | J| = 2 for ellipsoidal uncertainty of dimension 
20. Por testing the performance of maximizing the | linear optimization problems (6), we have 
created a single threaded C++ application that creates an affine function with coefficients randomly 
taken from the interval [—100,100], maximizes that function, randomly selects new coefficients, et 
cetera. Only the running time of the maximization step is measured. It can maximize 1.5 x 10® 
affine functions / : —)■ M per minute over a box, which is 10^® times slower than the MILP. It 

can maximize 3 x 10® affine functions g : —)■ M per minute over an ellipsoid, which is 32 times 

faster than the MIQCP. 

5.2 Illustrative small problems 

Consider the following toy optimization problem: 

(TOYl) min d 

s.t. d > max{x, X + C} + max{x, X — C} VCg[—1,1] 

d G M, X G M+. 


The optimal value of this problem is 1 (x 

= 0, d = 1). If we model this problem as an LP and then 

apply RO, we get the following model: 

(RC-R) min 

2/1 +2/2 



s.t. 

yi>x 

VC G [—1, 

1] 


2/1 > X + C 

VC G [-1, 

1] 


y2>x 

VC G [-1, 

1] 


y2>x-C 

yi,y 2 G M,x G M+. 

VC G [—1, 

1] 
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The optimal value of this model is 2 (x = 0, yi = y 2 = 1), so the RC-R is not exact. If we substitute 
yi = Vi + WiC, the model becomes: 

(AARC-R) min y 

s.t. y > vi + wiC + V 2 +W 2 C 

vi + wiC > X 

vi + wiC > x + C 

V2 + W 2 C > X 

V 2 + W 2 C, > X - C 
vi,V2,wi,W2,y G M,X G M+ 

The optimal value of this problem is 1 (xi = V 2 = ^,wi = ^,W 2 = —^,y = 1), which is the same 
as the optimal value of the original problem. So, in this case the AARC-R closes the gap. This is 
not always the case as the following example shows: 

(TOY 2) min d 

s.t. d > max{x, x -|- Ci + C2} + max{x, x -|- Ci — C2} 

-I- max{x, X - Cl + C2} + max{x, x - Ci - C2} VC G [- 1 , 1 ]^ 
d G M, X G M+. 

We obtain vtrue = 2, vrc-r = 8 and vaarc-R = 4. The AARC-R is an improvement over the 
RC-R, but is not exact. The QARC-R, which can be derived by reparameterizing the uncertainty 
region to [0,1]^, applying the work of Yamkoglu et al. (2012) to obtain an LMI description of 
the uncertainty region, and formulating the RC to a deterministic program by dualization, gives 
vqarc-R = 2.83. While this value is much closer to 2 than the value of the AARC-R, it is still 
inexact. 

5.3 Least absolute deviations regression with errors-in-variables 

An errors-in-variables regression model is a model y* = /3o + fdix* + s i (^£i N{0,a‘^) i.i.d.) where 

X* cannot be measured accurately. Only x* and yi with x* = (1 -|- Ci)xi are observed, where Q is 
an unknown measurement error. The least absolute deviations approach estimates l3o and l3i by 
minimizing Ya=i \yi - Po - Pix*[. 

n 

min Vmax{yj - Po - /3i(l + Ci)xi,Po + Pi{l + Ci)xi - yi}- 

P0,Pl • T 
1=1 

Because we do not know the values Ci; we can apply RO: 

n 

min max Vmax{yi - Pq - Pi{l + CpXi^Po -h Pi{l -h Qxi - yP}, 

/3o,Pi 

where we choose the uncertainty region Z to be ellipsoidal: ^ = {C G : IICII 2 — Note that 
this is the second special cases in Section 3.1.3 because it is the sum of absolute values, Q appears 
only in term i, and the uncertainty region is centrosymmetric around C = 0- h follows that the RC 
can be written as: 

n 

min Vmax{yj - Po - PiXi, Po -h PiXi - yi] + O ||^ix ||2 , (15) 

do,di 


VC e [-1,1] 
VC G [-1,1] 
VC G [-1,1] 
VC G [-1,1] 
VC G [-1,1] 
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which can be reformulated as a CQP. Even though there is this explicit and simple formulation of 
the exact RC, we also try the other approaches from Section 3 for a comparison. The RC-R is as 
follows: 


n 


(RC-R) 

min 

2 = 1 



s.t. 

Zi>yi- (do- (diXi -h Hj^iXij 

yiG I 



Zi> (do + (diXi -yi + Hj^iXij 

Vi G I 


which for comparison with formulation (15) can also be written as: 

n 

min Vmax{yj - /3o - PiXi, f3o + - yi} + O pixUi. 

/3o,/3i 

This is more pessimistic than constraint (15), since H-H^ > ||•|| 2 • For the RC-R, the worst case 
realization in the ellipsoidal set has all but one elements equal to 0. Therefore, it may seem that 
the RC-R is the RC of a linear constraint with interval uncertainty, but this is not the case. The 
AARC-R is given by: 


min d 



n 

s.t. d > ^ Vi + w^(( 

1 = 1 

VCGMV||C||2<f^ 


Vi + wlc>yi- (do- /3i(l -h Ci)xi 

VCGMV||C||2<f^ 

Vi G / 

Vi + wlC> /do + (diil + C,i)xi - yi 

VCGMV||C||2<f^ 

Vi G /, 


which after reformulation becomes: 

(AARC-R) min d 

n 

s.t. (j > n I |tt;| I 2 -|- Vj 

i=l 

Vi > Vi - l3o - PiXi -h n -h Wi\\2 yi G I 

Vi> /3o + PiXi -yi + 0,\\l3iXiei - Wi\\2 Vi G /, 

where e* is the i*^ unit vector. 

We have run these four models on 1, 000 cases. In each case, we fixed the parameters at /3o = 2, 
= 5 and = 1. A case consists of 15 observations of Xi, generated uniformly in [0,100] i.i.d., 
(( G : IICII 2 — ^ uniformly and Si ~ A^(0, a'^) i.i.d. Using these random draws, we have computed 
yi = (do + ,01 (1 -|- Ci)xi -|- Si- For the uncertainty region, we have picked U = 0.05. For each of this 
cases, we solved the exact formulation (15), AARC-R, QARC-R, RC-R, and the nominal problem 
using CVX and SDPT3. All solution times are very low, and hence, not mentioned. 

Next to comparing the usual way by means of vtrue, we also compare solutions by looking at 
how well /3o and /0i are estimated. In all 1,000 cases, the AARC-R and QARC-R give the same 
estimates of (do and ,0i. Hence, they are considered equal. Histograms, the mean, and statistics 
for f3o, (3i and vtme are shown in Figure 1. All models do well in estimating the parameters (do and 
,01, except the AARC-R and the QARC-R for /do- 

The objective values (averaged over the 1,000 cases) are listed in Table 1. The AARC-R and 
QARC-R objective values are 22.7% respectively 94.4% closer to the optimum than the RC-R ob¬ 
jective value. However, for the true value we have vtrueixRc-R) = 92.095 while vtrue{xAARC-R) = 


14 



Table 1: A comparison of different solution methods for the regression example averaged over 1,000 
datasets._ 


Method 

"^method 

'^true {^method) 

Nominal 

36.326 

91.994 

RC-R 

224.096 

92.095 

AARC-R 

194.091 

99.322 

QARC-R 

99.323 

99.322 

Exact 

91.973 

91.973 


vtrue{xQARC-R) = 99.322. The RC-R outperforms the AARC-R and the QARC-R in 98% of the 
cases, and even the nominal solution outperforms the RC-R in 99% of the cases. So, even though 
the AARC-R and QARC-R provide a much better bound on the optimal true value, their true 
value is not necessarily better. This shows two things: it may be very misleading to compare 
robust solutions by their objective values, and using a better approximation of the true objective 
function might not improve the solution at all. 

While we were able to solve the exact RC (15) directly, we have also looked at the efficiency of 
the cutting plane methods. This comparison is not based on the 1,000 cases mentioned earlier, but 
on new cases where we varied the dimension of C, (and consequently, the dimension of x). Each case 
was created in the same way as described before. Algorithms 1 and 2 can run very quickly because 
the worst case (( can be computed efficiently, as this is a simple case. It is therefore interesting to 
look at the number of iterations the algorithms take. We have generated 1,000 data sets in which 
we varied the number of observations. Algorithm 1 adds only 3 or 4 constraints, independent of the 
number of observations. This shows that the method is not only effective for polyhedral uncertainty 
regions. The number of iterations of Algorithm 2 is shown in Figure 2, where it seems to be a square 
root function of the number of observations. The regression model numitevi = a\/numobsi + e* 
gives a = 1.258 with a standard error of 0.007 and E? = 0.966. For 200 observations, the algorithm 
generates at most 18 constraints while full EORLC would result in a model with 2^^^ constraints. 
The algorithm needs more iterations than the vertex enumeration algorithm, but is still effective. 

5.4 Brachytherapy 

High dose rate brachytherapy (HDR-BT) is a form of radiation therapy where a highly radioactive 
sealed source is inserted into a tumor for short time periods via approximately fifteen till twenty 
catheters. When the catheter positions are fixed, a treatment plan specifies for how long the 
radioactive source has to stay at which position inside the catheters. A perfect treatment plan 
delivers a prescribed dose to the tumor while not delivering any dose to the surrounding organs at 
risk. Because this is physically impossible, the goal is to find a good trade-off between irradiating the 
tumor and saving the organs at risk. The quality of a plan can be measured by means of calculation 
points. These are artificial points where the received dose can be computed and compared to a 
prescribed lower and upper bound, that are distributed inside and around the tumor and organs 
at risk. The dose in a specific calculation point i can be computed as the sum of the individual 
contributions from every catheter k, which in turn is the sum of the dwell times of the individual 
dwell positions inside catheter k multiplied by given dose rate vectors dik- 

k&K 
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Figure 1: Comparison of the nominal solution, the RC-R solution, the AARC-R/QARC-R solution 
and the exact solution for the regression model of Section 5.3. 




Number of data points 


Figure 2: The number of iterations needed by Algorithm 2 to find an optimal solution versus the 
number of observations in the regression model of Section 5.3. 
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If calculation point i fails to satisfy the lower bound Lj or the upper bound Ui, it contributes a 
linear penalty of ai or /3i (respectively) per unit of violation to the objective function. This results 
in the following optimization problem: 

d]ktk),(3i{Y d]ktk - Ui)}. 

This convex piecewise linear objective function is commonly used for treatment planning (Alterovitz 
et al. (2006); Karabis et al. (2009); Lessard and Pouliot (2001)). 

The parameters dik are computed based on the catheter positions. These positions cannot be 
measured accurately, hence the data in the optimization problem is uncertain. We assume that the 
true catheter position is within some cone around the measured position, which is justified because 
one side of the catheter is fixed at a known position. We replace the cone with a polyhedral cone 
with a 10-sided base to end up with a polyhedral set, and we make the simplification that the 
vector dik inside the cone is a convex combination of the vectors dik at the sides of the cone. Thus, 
we only need to know the vectors dik at the catheter positions corresponding to the 10 edges of the 
cone connecting the apex to the base. Using these vectors as the columns of a matrix B, we can 
write dik = BikCk, where Cfc G = {C £ : e^C = 1,C > 0}, the standard simplex 

in RI'^'I), and |5| is the number of sides of the base of the polyhedral cone. This gives the following 
RO problem: 

min V 

s.t. Y. maxjO, ai{Li — E {BikCkf tk), 

i k&K 

A( E iBikCkYtk - Ui)} <v va G aI^I-^A: G K) 

k&K 

tfe > 0 yk e K. 

The number of calculation points is usually in the order of magnitude of 5,000. We have data 
of one patient that has been exported from treatment planning software with a large number of 
calculation points. For the purpose of this paper, we have reduced the number of calculation points 
to 40 by taking a random subset, because otherwise both the maximization step in Algorithm 2 and 
the AARC-R are intractable. The model does not lose its value from this reduction, because the 
objective could be a weighted average between the nominal objective based on all calculation point, 
and the robust objective based on a small fraction of calculation points that is well distributed 
inside the tumor. 

The RC-R is based on the original constraint with the sum of 40 maxima of 3 functions. We have 
also solved the RC-R and AARC-R after splitting up the sum, as outlined in the last paragraph of 
Section 3.2. If we split the sum in groups of size 4, then the number of sums reduces to 10 while 
each term of sum is the maximum of 3^ functions. 

The results for the different methods are listed in Table 2. The first observation is that the 
nominal solution has a true value that is five times the optimal value, so the nominal solution is 
nonrobust. The AARC-R has value ~ 143, which is much closer to the true value (~ 137) than 
to the RC-R value (~ 229), so the AARC-R almost closes the gap. Algorithm 2 is outperformed 
by Algorithm 1. Its long running time when the optimization problem for determining Vtme is 
not solved to optimality is not only due to the larger number of iterations, but more importantly, 
to the gradually increasing running time of solving LP in each iteration. This step becomes so 
memory consuming, that the algorithm cannot finish with CPLEX on 32 bit hardware. Combining 
both algorithms reduces the number of iterations, as more work per iteration is done, but also 
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Table 2: A comparison of different solution methods for the brachytherapy example. 


Method 

Iterations 

Sol. time (s) 

'^method 

"^true {^method) 

Nominal 

— 

0 

11.033 

711.362 

RC-R 

— 

1 

229.986 

157.277 

AARC-R 

- 

314 

143.104 

140.827 

RC-R (sum of 8 maxima of 243 functions) 

— 

51 

186.219 

157.568 

RC-R (sum of 10 maxima of 81 functions) 

— 

17 

203.100 

157.030 

RC-R (sum of 20 maxima 9 functions) 

— 

2 

218.894 

158.438 

AARC-R (sum of 10 maxima of 81 functions 

— 

136214 

141.481 

140.114 

AARC-R (sum of 20 maxima 9 functions) 

- 

2774 

142.447 

139.957 

Algorithm 1 

25 

1,173 

137.628 

137.676 

Algorithm 2 

36 

1,065 

137.629 

137.662 

Combined 

20 

894 

137.628 

137.650 

Algorithm 1* 

53 

234 

137.627 

137.668 

Algorithm 2* 

> 400 

> 50,000 

- 

- 

Combined* 

41 

353 

137.628 

137.676 


e = 0.05 

* The optimization problem for determining vtrue is stopped as soon as the upper bound is at least 0.05 higher than the lower 
bound 


increases the solution time in case the optimization problem for determining vtme is not solved to 
optimality. Splitting up the sum does not give any benefit in terms of true value, though it gives a 
lower objective value when solving the RC-R. 


5.5 Inventory planning 

We consider a single item inventory model where backlogging is allowed to compare the AARC- 
R with the exact RC. At the beginning of each period, an order can be placed that is delivered 
instantly. At the end of each period, the holding and backlogging costs are Ch and Cb per unit, 
respectively. The objective is to minimize the costs: 

/m t 

min max ^max{c/i[xo + E Qi - di],Cb[xo + 

\t=l i=l 



where xq is the starting inventory, qi is the order quantity at time i, and di is the uncertain demand 
at time i. In all formulations we allow to depend affinely on the demand in periods 1 up to i — 1. 

If the uncertainty region is a box, the AARC-R turns out to be exact for the intervals we 
have tried. This is in accordance with the numbers reported by Ben-Tal et al. (2005). We found 
that the AARC-R is no longer exact if the uncertainty region is an ellipsoid. Because demand is 
nonnegative, the ellipsoid is intersected with the nonnegative orthant: 


Z = {de 


d 



< II}. 


In our numerical study we have looked at 12 time periods, with parameters H = 10, d = 5, Ch = 1, 
and Cb = 2. 
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Each determination of vtrue in Algorithm 2 takes up to 2 minutes if continued to optimality 
using CPLEX, making it the most time consuming step in the algorithm. Because the dimension 
of the uncertain demand vector is 12, a global solver might be faster. We have tried LGO 1.0, 
whose accuracy can be adjusted with the parameters “maximal number of function evaluations" 
and “maximal number of stalled evaluations", which initially both are 16,000. Every time the upper 
bound found by LGO is less than 0.1 larger than the lower bound, the parameters are increased 
by 25% until they exceed 1,000,000. This reflects the idea that it is still easy to find a violated 
constraint when the algorithm starts, but gradually becomes more difficult as the quality of the 
solution increases. Still, GPLEX finds better solutions in less time. 

The RG-R is based on the original constraint with the sum of 12 maxima of 2 functions. We 
have also solved the RG-R and AARG-R after splitting up the sum, as outlined in the last paragraph 
of Section 3.2. In this example the uncertainty enters the model through time, so the most natural 
way of splitting up the summation is in consecutive time periods. We split the problem into two 
groups: 


min yi + y 2 

6 t t 


s.t. 

E 

max{ch[xo + '^qi- di], Cb[xo + '^qi-di] < yi 

^d^Z 


t=i 

i=l 

i=l 



12 

t 

t 



E 

t=7 

max{ch[a:o + E “ '^*1 ’ ^b[xo + E 

i=l i=l 

yd£Z 


where again we allow qi to depend affinely on the demand in periods 1 up to i — 1. This reformulation 
introduces more constraints with the sum of maxima, but each sum contains less terms, hopefully 
resulting in a shorter solution time. The worst case d may differ for the constraints that set yi and 
y 2 , hence this reformulation is not exact. 

The results are listed in Table 3. Again note that the order policies {qi) are adjustable in 
all formulations, including the RG-R, so the differences in this table are caused only by different 
reformulations of the sum of maxima. The first observation is that the AARG-R gives a very small 
gain over the RG-R, but still has almost twice the optimal value. So, making the analysis variables 
adjustable does not significantly improve the solution. Algorithm 2 is the fastest cutting plane 
method, requiring approximately the same number of iterations when combined with Algorithm 1. 
The sum splitting method significantly reduces the computation time at the cost of nonoptimality. 
It performs much better than the AARG-R, both in optimal value and in true value. Using the 
AARG-R on the splitted sums gives a very small gain over the RG-R on the splitted sums, just as 
for the full problem, so it is not listed in the table. The large true value of the nominal solution 
comes from the fact that the order sizes are fixed in advance and are not adjusted to observed 
demand. 

6 Conclusions 

Because RO is applied constraint-wise, it is very important how constraints are formulated. In 
this paper we list several approaches to an inequality constraint containing the sum of maxima of 
linear functions. The RG-R, often used in the literature, is the most pessimistic approach. It is 
obtained by first reformulating the deterministic constraint into linear constraints using analysis 
variables, and then applying RO. Its pessimism can be reduced by replacing the analysis variables 
with linear decision rules before applying RO, which gives the AARG-R. The AARG-R seems to 
work well for the practical problems we analyzed with polyhedral uncertainty regions, but we have 
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Table 3: A comparison of different solution methods for the inventory example. 


Method 


Iterations 

Sol. time (min) 

^method 

'^truei^method) 

Nominal 


— 

0 

0.000 

509.903 

RC-R 


— 

0 

120.000 

94.641 

AARC-R 


- 

0 

120.000 

93.315 

EORLC 


— 

10 

48.750 

48.750 

RC-R (sum of 2 

maxima of 64 functions) 

— 

0 

68.613 

54.162 

RC-R (sum of 3 

maxima of 16 functions) 

— 

0 

83.631 

58.140 

RC-R (sum of 4 

maxima of 8 functions) 

— 

0 

94.456 

66.637 

RC-R (sum of 6 

maxima of 4 functions) 

- 

0 

107.627 

77.957 

Algorithm 1 


123 

76 

48.749 

48.796 

Algorithm 2 


77 

36 

48.752 

48.802 

Combined 


79 

45 

48.755 

48.798 

Algorithm 1* 


309 

180 

48.752 

48.827 

Algorithm 2* 


112 

9 

48.753 

48.790 

Combined * 


103 

29 

48.755 

48.825 


e = 0.1 


* The optimization problem for determining vtrue is stopped as soon as the upper bound is at least 0.1 higher than the lower 
bound 


constructed an example with polyhedral uncertainty where the value of the AARC-R is 100% higher 
than the true value. Nonlinear decision rules may give better results, but are computationally more 
challenging. The conservatism of the approximations can be reduced by combining several max 
expressions before reformulating. Especially for ellipsoidal uncertainty this method gives much 
better solutions at the cost of a slightly higher solution time. 

In many cases it is not necessary to use an approximation because an exact reformulation can be 
practically solved. We identify four special cases in which an exact reformulation is often tractable. 
For the general case we give two exact general methods: vertex enumeration and EORLC. Both 
methods may result in very large optimization models, but cutting plane methods can be used to 
handle this. Vertex enumeration adds a set of constraints for every vertex of the uncertainty region, 
so this method is preferred if the uncertainty region has a low number of vertices. Surprisingly, 
its cutting plane version is also capable of solving problems where the uncertainty region has an 
infinite number of extreme points efficiently. EORLC is preferred if the number of terms with a 
maximum function is low. If it is not clear in advance which cutting plane method is faster, both 
methods have to be tried because our numerical examples do not show a clear preference. Both 
methods can be combined, but we have not found a situation in which it is beneficial to do so. 
EORLC can be combined with approximations, which is a new idea that can be used for large scale 
problems which shows promising numerical results. 

The RC-R is often used in the literature while less conservative approaches could have been 
applied, mostly without explicitly mentioning that their approach is conservative. In the paper 
by Kropat and Weber (2008), the exact method EORLC would have increased the number of 
constraints by only a factor four while reducing the number of variables with almost a factor 2 
and not changing the structure of the problem. The same authors applied vertex enumeration to 
an ellipsoidal constraint with polyhedral uncertainty (Ozmen et ah, 2011), which gives an exact 
reformulation. Bertsimas and Thiele (2006), Wei et al. (2009) and Alem and Morabito (2012) apply 
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the RC-R to a problem with polyhedral uncertainty. Our results, and also the numbers reported 
by Ben-Tal et al. (2005), show that the AARC-R often gives less conservative solutions while it has 
the advantage that the problem remains linear. For inventory problems in general, when the order 
quantities are made adjustable, then oftenly also the analysis variables are made adjustable. We 
show that the latter is not beneficial for ellipsoidal uncertainty, and that both the RC-R and the 
AARC-R give very bad solutions. For small planning horizons, an exact method has to be used, 
while for larger horizons splitting the sum in small groups and applying EORLC on the groups 
significantly improves the solution. Ng et al. (2010) solve a lot allocation problem with ellipsoidal 
uncertainty. Because the problem is computationally challenging, they solve a problem equivalent 
to our RC-R using Benders decomposition. Even though their problem is so challenging that even 
the simple RC-R cannot be solved within a day and their speed-ups are necessary to solve the 
problem efficiently, it is still interesting to know the conservatism of their approach. We have been 
able to get a suboptimal solution with cutting planes based on vertex enumeration, where both the 
minimization and the maximization step were stopped before optimality. The solution we got after 
four hours has a true value of 26.8, whereas the RC-R (which we tried to solve as an MILP) has a 
value between 28.7 and 47.2. So the objective value of the solution proposed by Ng et al. (2010) is 
at least 7-76% too pessimistic. 

From our numerical examples it becomes clear that the RC-R is not necessarily better than the 
nominal problem. Neither of the two optimizes the true problem, so it cannot be determined a 
priori which one has a better true value. The same holds for the AARC-R and the RC-R: If the 
AARC-R gives a much lower value then at least it provides a guarantee on the worst case, but 
the RC-R may still outperform the AARC-R. When using an approximation, it is therefore crucial 
to measure its quality. This can be accomplished by comparing the true value of the solution of 
the approximation with the value of an exact formulation. If the problem is too large to be solved 
exactly, the comparison may be based on a smaller instance with similar structure. 
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A Derivation of AARC-R using Fenchel’s duality 

In this appendix we apply Fenchel’s duality to robust constraints, a technique introduced in RO 
by Ben-Tal et al. (2014). First we will briefly mention the general theory, then we will apply it to 
constraints of general form, and finally we will apply the general results to constraint (1) and show 
that the result is the same as the AARC-R. 

A.l Fenchel’s duality theorem 

We start with some definitions that are necessary to formulate Fenchel’s theorem: 

Definition A function (p is proper convex if it is convex, its codomain is MU {oo}, and (p{x) < oo 
for at least one x. 

Definition A function ■0 is proper concave if it is concave, its codomain is MU {—oo}, and 0(x) > 
—oo for at least one x. 
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Theorem A.l (Fenchel’s duality (Rockafellar, 1970, p. 327)) Let c) be a proper eonvex funetion 
on M”, let be a proper eoneave funetion on M”, and let either of the following eonditions be 
satisifed: 

• ri(dom 4>) Cl ri(dom if) 

• 4> and xjj are elosed, and ri(dom (f>*) n ri(dom 'if^) / 0, 

where ri is the relative interior, doni is the effeetive domain (dom f) = {x ■. 4>{x) < oo} ), and (f* 
and V'* are the eonvex and eoneave eonjugate of (f and fj, respeetively. That is, 

(t)*{s) = supjs^x — 4>{x)} 

X 

Then the following equality holds 

inf{(/)(x) - ■0(x)} = sup{V^*(s) - ())*{s)}. 

X g 

A.1.1 Fenchel’s duality applied to a robust constraint of general form. 

We focus on the general robust constraint 

g{C,x)<d yfeZ, (16) 

where (/ is a proper concave function of f for any fixed value of x, and the condition for Fenchel’s 
duality is satisfied with respect to the first argument for any fixed value of x. Because values of g 
are not of interest when f ^ Z, we assume that g{(), x) = —oo for all f ^ Z. We also assume that 
Z \s Si compact set so that this constraint is equivalent to: 

max{5r(C, x) - < 5 ^( 0 } < d, (17) 

with 5z the indicator function {5z{C} = 0 if ( £ Z, and oo otherwise). We can rewrite the left-hand 
side by applying Fenchel’s duality: 

min{(i^(s) — g^{s, x)} < d, 

which holds if and only if there exists some s £ such that: 

d*z{s) - g*{s,x) < d. (18) 

Note that this constraint is convex in s. Because every step is ‘if and only if’, constraint (16) is 
equivalent to constraint (18). 

A. 1.2 g is the sum of other functions. 

If g can be written as the sum of several other functions, it might be impossible or very difficult to 
find a closed form solution for its conjugate function. Suppose we have a constraint of the form: 

Y.9i{C,x)<d yc^z, (19) 

i&I 

which is constraint (16) with g = gi- If we want to formulate an equivalent constraint using 
Fenchel’s duality, we need the concave conjugate of g. Under some mild assumptions on gi, it turns 
out to be sufficient to have closed form solutions for the conjugates of gi. The following lemma 
appears to be very useful (Rockafellar, 1970, p. 145): 
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Lemma A.2 Let (i G I) be proper eoneave funetions on M"'. //Pljg/ ri(dom 'ipi) / 0 then 


{^'4^i)*{s)= sup 
iel i&I 


and the supremum is attained. 

Applying this lemma, we can rewrite constraint (19) as: 

- max a:)} < d, 

i&l 


which is valid if and only if there exists s* G M^(i G I) such that 

dziYl '®*) “ x) < d. (20) 

iel i£l 

So, if all conjugates have closed form solutions, the reformulated constraint also has a closed form. 
If all conjugates do not have a closed form, this reformulation is still useful because it allows 
computing each term separately. We will show this for piecewise linear convex functions later in 
this section. It should again be noted that this constraint is convex in Sj. 


A.1.3 / is not convex or g is not concave. 

Let us investigate the consequences to the reformulation if g is not concave. We assume g is finite on 
some non-empty set and oo elsewhere, so that its conjugate is non-trivial. Fenchel’s inequality 
(Rockafellar, 1970, p. 105) states that: 

dziC) + d*z{s) = SziC) + sup{C'% - dziC')} > dziC) + - 5z{C) = s 

C'&z 

g{C,x) + g*{s,x) = g{C,x) -h mf {C'% - g{C',x)} < g{C,x) -h C^s - g{C,x) = Cs, 

hence: 


and consequently: 


dz{C,) + d*z(s) > g{C,x) + g4s,x), 


g{C, x) - dz{C) < d*z{s) - g^{s, x). 

This implies that if constraint (18) is satisfied, then so is constraint (17), but the reverse implication 
is not necessarily true. Hence, constraint (18) is a conservative reformulation of constraint (17). 
Also, Lemma A.2 no longer holds with equality. We can rewrite it as an inequality: 

= inf{s^t -J^^i(t)} 
i£l i&I 

= inf{^sJt-V’i(t)}, 

i£l 
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for any si ,Sm £ for which = s. So in particular: 

i&I Y.iei^i=^ *e/ 

> sup {^inf{s-t - V'i(t)}} 

is/ 

= sup {Xl(V’i)*(s*)}- 

*e/ 

This implies that if constraint (20) is satisfied, then so is constraint (19), but the reverse implication 
is not necessarily true. Hence, constraint (20) can be seen as a conservative reformulation of 
constraint (19). 


A.1.4 g is the sum of pointwise maxima of linear functions. 

Let us first derive the conjugates of some functions before we arrive at the theorem: 

5z{s) = sup {s\ - fe(C)} = supjs^C} = max{sY}, 

iez Q&z 


and 


(max{£i,(C,x)})*(sj,x) 

j&J 


inf {s-C - max{fjj(C,x)}} 


C6.Z 


jeJ 


inUminis] C - ^ij{C,x)}} 
C,ez jeJ 


= mini inf jsJC 
jeJ 




Theorem A.3 Applying Fenchel’s duality to: 

max 
C6.E 

gives a formulation that is equivalent to the AARC-R. 


Vmax{£ij(C,x)} < d, 

- r 

zei 


( 21 ) 


Proof Constraint (21) is equivalent to constraint (19) with gi{C,x) = maxjgj{fjj((^, x)} for C, in 
Z and gi{f,x) = —oo otherwise. For for any fixed x, gi is not concave in f so we will end up 
with a conservative instead of an equivalent reformulation. If we fill in the conjugate functions in 
constraint (20), the following conservative reformulation is obtained: 

If we model the second terms as ^an write this as: 

'^slC-'^Zi^d \lf^Z 

i£l i£l 

Zi < s]( — iijic, x) e z yi e I Vj G j, 

and by rearranging the terms in each constraint we obtain: 

^[(—^i) + SjC] < <31 yf £ z 

iei 

{~zi) + Sjc ^ ^ijiCix) yf £Z yi £ I yj £ J, 

which is the same as the AARC-R. | 
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B Derivation of AARC-R by reformulating the nonrobust con¬ 
straint 


In this appendix we give a different derivation of the AARC-R of constraint (1) when both the 
biaffine functions and the uncertainty region are separable in the following way: 

^max{^ < d {k G K), (22) 

i&I keK 

and Zk is the convex hull of different scenarios {s G S). An example where this constraint 
is commonly used, is HDR brachytherapy optimization (Alterovitz et al. (2006); Karabis et al. 
(2009); Lessard and Pouliot (2001)). If the summation over k were outside the max expression, 
then an analysis variable could be used for every k without introducing any conservatism. Vertex 
enumeration can then be done on every Zk separately. For problems not affected by uncertainty, 
we show that indeed there is an equivalent formulation where the summation over k is outside the 
max expression. Then we show equivalence to the AARC-R if there is uncertainty. First, we prove 
the following equality for fixed x and (k- 

Lemma B.l 


H max{ ^ 
iGl k^K 


mm 


k^K 


Vijk —0 


EE 

k&K iel 


max{yijk + iijkiCk, x)}. 

j&J 


(23) 


Proof Note that: 


EE max{yijk -h iijkiCk, a^)} > EE yij(i)k “1“ ^ij{i)k{Ck^ EE hj(i)k{Ck,x) yj{i) G J, 

k&K i&I k&K i£l i&I k&K 


for any j(i) in J, so in particular: 

EE max{yijk -h iijkiCk, a:)} > E max{ E ^ijkiCk^ a:)}, 

k&K i£l i&I k&K 

for any y, so in particular for the minimum. Hence, the right hand side of (23) is at least as large 
as the left hand side. On the other hand, given a feasible point for the left hand side of (23), we 
can always construct a feasible point for the right hand side with equal value by taking the same 
X, and yijk — 'Y^k'&K ^ijk'iCk'iX) (-ijkiCkj x)'. 

E + ^ijk{Ck,x)} = E E E ^ijk'{Ck',x)} = ^ max{ ^ iijk'iCk', x)} ■ 

k£K i£l i&I k&K k'&K i&I k'&K 

I 

Lemma B.2 A conservative reformulation of constraint (22) is given by: 

V V max{yijk -h iijkiCk, x)} <d VCfc G Zk {k G K) 

E yijk = 0 \/iGl \/jGJ. 

k&K 
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Proof We use the proof of Lemma B.l. The hrst part of the proof still holds, so the equality in 
constraint (23) becomes a The construction of the feasible point in the second part of the 

proof depends on the value iijkiCk, x). This value does not exist in the robust constraint (22), since 
there is no single (k- Note that the left hand sides of (22) and (23) are the same, so by replacing the 
left hand side of (22) with the right hand side of (23), the conservative reformulation is obtained. 

I 


The uncertainty is separable per k in the first constraint of the conservative reformulation. Hence, 
analysis variables per k can replace each term without changing the solution. Vertex enumeration 
can then be done for every Zf^ separately: 


^ Zk<d 
k&K 

^k — ^ ' ^iks 
i&I 

Wiks Z Vijk T ^ijki^Ckt x) 
Vijk — 0 

k&K 


(24) 


Vs G 5 

Vi G / Vs G 5 VA: G iV 
Vi G / Vj G J. 


It remains to show that this formulation is equivalent to the AARC-R of constraint (22). 

Theorem B.3 The conservative reformulation in Lemma B.2 is the AARC-R in case of scenario 
generated uncertainty. 


Proof The AARC-R of constraint (22) is given by: 


^ Ui ^ w]kCk < d 

iei \ k&K / 

Xi + 'y ) '^ikCk ^ y ' ^ijiCk-i 

k&K k&K 

Constraint (25) can be reformulated as: 

wlkCk > yijk + iij{Ck,x) 

y ' Uijk — Xi 
k&K 


'iCk(^Zk {k£K) 

x) yiel Vj G J VCfc G Zk {k G K). (25) 


Vi G / Vj G J VCfc ^Zk ykeK 
yiel Vj G J. 


We have assumed scenario generated uncertainty, so our uncertainty region is Zk = tlie 

standard simplex in Suppose an optimal solution has Vi 0, then an optimal solution with 
Vi = 0 can be obtained by increasing all elements of wik for a single random k with Vi because the 
elements of fk sum to 1. Hence we fix Uj = 0. The AARC-R can then be formulated as: 


Zk<d 


k&K 


Zk>Yl ^ikCk 


iei 


wlkCk > Vijk +^ij{Ck,x) 
Vijk = 0 

k£K 


VCfcG^fc {k£K) 

yiel yj eJ VCfc eZk yk e K 
Vi G I Vj G J. 


Let Wiks denote the component of Wik- Equivalence to (24) now follows from vertex enumeration. 
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C Derivation of the QARC-R for an ellipsoidal uncertainty region 


In this appendix the SDP reformulation of the QARC-R of (1) is derived for an ellipsoidal uncer¬ 
tainty region. For simplicity, we assume iij{C,x) and i are bilinear functions in the parameters. 
Therefore, they can be expressed as £(C,x) = Lx and £ij((,x) = (^Lijx, respectively, for some 
matrices L and Lij. The QARC-R is given by: 

(QARC-R) C'^a: + ^(ui + u;-C + C'tFiC) <d VC £ : IICII 2 < 

iei 

Vi + w]C + CWiC>CLijX VCgM^ : IICII2 Vi G / Vj G J. 


By application of the 5-lemma (Ben-Tal et ah, 2009a, Lemma 6.5.3), this can be reformulated as 
the following LMIs: 




i&I \2 


Wi Vi 


Wi 


^^Lx —d ^ 

\ {wi - Lijx) 


2\'^i 

I {wi - Lijx) Vi 


^0 -Q?j 

- 0 -02 


Ao > 0,Aij > 0. 

The QARC-R has |/|| J| -M LMIs of size |L -h 1|. 


Vi G / Vj G J 
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